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The complete system of the B-spline solutions for the Dirac equation with the parity- 
nonconserving (PNC) weak interaction effective potential is obtained. This system can be used 
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matrix elements with relative accuracy 10 _r . 

PACS number(s): 31.30.Jv, 31.10.+Z 



I. INTRODUCTION 

Consideration of the parity nonconservation (PNC) effects in atoms provides an important verification of the 
Standard Electroweak Model in the low energy sector. The most important information provides the neutral Cs 
atom, because the experimental and the theoretical accuracy is the best in this case. The analysis performed in [1] 
and [2] has indicated a deviation of the measured weak charge value Q w from that predicted by the Standard Model 
by 2.5 standard deviations a. Later this value was diminishing and growing up again in a series of works [3]- [7]. 
It was understood that the radiative corrections to the PNC amplitude play an important role. This question was 
investigated in [8]- [12] where both electron self-energy (SE) and vacuum polarization (VP) corrections were evaluated 
by different methods. However, in these calculations aZ expansion for SE (a is the fine-structure constant, Z is the 
charge of the nucleus), Uehling approximation for VP or other approximations were employed. To our mind the direct 
accurate calculations of SE and VP without the use of any approximations are desirable. For this purpose we propose 
in this paper the accurate spline solutions of the Dirac equation with the weak-interaction PNC potential. 

Another possible area of application of PNC spline solutions could be the spectra of the Highly- Charged Ions (HCI). 
Though no experimental results are available here up to now, the analysis of the PNC amplitude in the two-electron 
ions can also provide an important test of the Standard Model [13]- [19]. The calculation of correlation effects in 
two-electron ions is much easier than in neutral atoms, but QED corrections are more important. 

The rigorous direct way to the calculation of SE and VP corrections to the PNC amplitude (see [20]) is rather 
cumbersome. In order to avoid the difficulties and receive the accurate results we propose to use another way: 
accurate solutions of the Dirac equation with the weak-interaction PNC potential. Having got the solutions of the 
Dirac equation with this potential we can obtain self-energy corrections, using the standard methods [21] - [23] based 
on the spline approximation [24], [25]. In this work we solve the first part of the problem: the obtaining of the 
complete basis set of accurate solutions of the Dirac equation with the weak PNC potential. The main difficulty 
is the extreme smallness of the weak interactions in the atomic scale. To overcome this difficulty we change the 
weak interaction potential by scaling it in such a way that it becomes significant but still small enough to apply 
the perturbation theory. Then we can obtain the necessary PNC matrix elements by simple rescaling. Actually our 
numerical procedure appear to be so accurate that the scaling will become necessary only in few cases, for the most 
singular operators and the enhancement of the weak interaction potential will not exceed 10 times. 

II. DIRAC EQUATION WITH PNC WEAK INTERACTION 

We employ atomic units m e = e = H = 1 throughout the paper. Consider stationary Dirac equation with the 
weak-interaction PNC potential 

[cap + V(r) +[3c 2 + V w (r)] ^{r)=E^{r) (1) 

where a, (3 are Dirac matrices, p, E are the electron momentum and energy, c ~ 137.03599976 is the speed of light, 
V(r) is the Coulomb potential of atomic nucleus (pointlike or extended) and V w is the weak-interaction PNC potential. 
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This potential (see, for example, [26]) looks like 
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G F 
2V2 



Qw Pnuc 

(r) 



(2) 



75 



I 

1 



Q w = - N + Z (1 - 4 sin 2 6 W ), 
Pnuc{r) = po [1 + exp[(r - c) / a]' 
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where 6 W is the Weinberg angle, sin 2 9 W w 0.2230 [27], iV is the number of neutrons in the nucleus, Z is the number of 
protons, Gf is the Fermi constant and p n uc{r) is the Fermi distribution for nuclear charge density. The eigcnfunctions 
of the Dirac equation with PNC potential have no definite parity, therefore we assume 
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where ip ( r ) and % ( f) are the upper and lower components of Dirac bispinor, ( r ), g 2 ^ ( r ), ( r ) and / 2 j( ( r ) 
are the radial functions, fljiM and f^f M are spherical spinors with opposite parity. The equations for the components 

<P, X are 



(E - V{r) - c 2 ) cp(r) - c(ap) X (r) - W(r) X (r) 



Inserting (4) in (5), (6) and separating out the angular variables we get 
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This system of equations for the radial functions was first derived in [28] . However in [28] only one certain solution 
of Eq (7) was obtained. Our finite goal requires the knowledge of the complete system of the solutions of Eq (7). For 
this purpose we use the B-spline approach [24], [25]. 



III. PARITY-NONCONSERVING B-SPLINE APPROACH 



Using the Galerkin method (see [29] and [24]) we express the system of equations (7) in terms of the least action 
principle. We use MIT bag- model boundary conditions [30] , employed also in [24] . In case of PNC functions these 
conditions look like 

g 1 (R) = f 1 (R) fff 1 (0) = 

g*(R) = f 2 (R) \ 3 2 (0)^0. U 

Here R is the size of a spherical box, where an atom or an ion is enclosed. We choose action functional S so that 
condition 5S = leads to Eqs (7) and (9). Then we present the functions g 1 , g 2 , f 1 and f 2 in the form of the linear 
combination of B-splines. The condition 5S — reduces to the system of 4n x 4n symmetric generalized eigenvalue 
equations. For constructing the B-spline system we use the grid r, = pj/Z, i = 0, . . . , N int , that was described 
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in [25]. In our calculations we use the number of the grid intervals Ni nt — 200, the order of B-splines k — 9 and the 
parameter, which defines the density of knots near zero (see [25]), h se t = 0.0167. 

In this paper we investigate the hydrogenlike ion 238 U 91+ . The box size was taken to be R « 1.0 a.u. The parameters 
of the Fermi distribution arc c = 7.136 fm and a = 2.3/4 In 3 [31]. As a first test of the accuracy of our B-splinc 
approach, we set V(r) = —Z/r, N w = and compare the energy values from B-spline approach with low-lying energy 
levels given by the Sommerfeld formula. The results are presented in the Table I in a.u. Our accuracy of calculation 
of the energy 10~ 14 is the limit of accuracy for the calculation with Fortran double precision type of variables. 

Another check of the spline accuracy is given in Table II, where the matrix elements of the operators r m (m — 
2, 1, 0, —1, —2) arc compared on the 2s\j 2 and 2p 1 / 2 spline and exact (point-like nucleus) wave functions without weak 
interaction. The results in Tables I and II show that all the matrix elements < 2s 1 / 2 |r m |2s 1 / 2 > and 
< 2p 1 / 2 |r m |2p 1 / 2 > apart from the case m = —2 can be evaluated with the same relative accuracy 10~ 14 — 10~ 15 as 
the low-lying energies. Only in case of m = — 2 the accuracy is of order 10~ 7 — 10~ 8 because of the singularity of r~ 2 
operator. 

IV. SCALING AND THE NUMERICAL TESTS 

As the next step we increase the PNC weak interaction by scaling V w = N w r )zW(r) 1 where N w is the scaling 
factor, and investigate the energy spectrum by changing N w . Here we use the Fermi distribution for the nuclear 
charge density. 

The scaling problem looks as follows: to find a scaling parameter N w for which 1) the contribution of weak interaction 
to the matrix elements of our interest is larger than numerical inaccuracy, 2) weak interaction is still small enough 
for using the perturbation theory. 

The Table III presents the PNC weak interaction contribution to the energy of the states for the different values of 
N w . In case N w ^ the states do not have a certain parity and we use the notations ns, np formally classifying the 
different states by their origin. The analysis of the results of calculations shows that the dependence E = E(N W ) is 
approximately parabolic. One can also obtain this result from the perturbation theory (PT): 

|n S ~) = |n S ) + N W £ ^f^ns) ^ + 1 ( ^ )2 ^ ^ ^V w \rnp}(^p\V w \ns} (iQ) 
(ns\H\n~s) - E ns = - (N w f £ (11 ) 

l^ns ^mp 

Here the evident equalities (ms|K|ns) = and (mp\H(r)\ns) = were used. The right-hand side of Eq (11) is also 
presented in Table III for N w = 1 . The comparison shows that we can obtain the PNC weak interaction contribution 
to the energy with accuracy 0.03 % if we use scaling parameter N w = 10 2 -j- 10 3 . Here by the numerical inaccuracy we 
undestand the B-splinc method inaccuracy of calculation of nonperturbative energy value presented in the last column 
of Tabic III. This leads to the poor accuracy of the weak-interaction PNC contributions to the energy, obtained with 
the low values of scaling parameter N w = 1 10. For N w > 10 4 higher orders of perturbation theory become to be 
not negligible. 

The second-order PNC weak interaction correction to the energy is of no physical interest and was investigated 
here only for the demonstration of the possibilities of the scaling procedure. This procedure helps to evaluate the 
second-order PNC weak interaction correction even though this correction is beyond the level of accuracy of the spline 
approximation. 

We should emphasize that the matrix elements of our interest are the matrix elements with linear dependence on 
weak interaction. SE and VP corrections to the PNC amplitude are the objects of this type. 

For testing our system of PNC B-spline functions we calculated the set of the matrix elements (2s|r m |2s) with the 
first non-zero order of perturbation expansion, linear in V w . The following tests were employed: 



(og| r m |2p) = V V 1 [lifA^fll | \npr<ns\r ■ \zp/ 

~f E 2 p — E n s ~f E 2 p — E n p 



[<2~\r m \2P) = Y I V i^™gJ^ (13) 

' E 2n — E ns J Els — Enn 
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With these tests we check the completeness of the PNC B-spline spectrum (Eq (12)) and the absence of the higher 
order PNC contributions after the introduction of the scaling parameter N w i.e. the possibility of rcscaling (Eq (13)). 
The results are given in Table IV. The behaviour of the matrix elements with m = 2, 1 is similar. The most suitable 
scaling parameters for the calculation of these matrix elements are N w = 1 (the absence of the scaling) with relative 
accuracy 1 • 1CP 8 and N w = 10 with relative accuracy 5 • 1CP 8 . The test of completeness of the spectrum for these 
operators gives the relative accuracy 1 • 10~ 12 . The most interesting matrix elements are the ones with m = —1,-2 
because of their singularity. For N w = 1 the relative inaccuracy for the calculation of both matrix elements is about 
1 • 10~ 6 , for N w — 10 it is 1 • 10~ 7 . The checking of completeness of the basis for this operators give the relative 
accuracy 1 • 10~ 7 in the worst case. Summarizing these tests we conclude that the best scaling parameter for the 
calculation of matrix element with linear dependence on weak interaction is N w = 10. The accuracy of this calculation 
is about 1 • 10~ 7 . 

Finally, we want to illustrate our method by calculation of VP corrections in the Uehling approximation to the 
matrix element of PNC weak interaction on \2s l / 2 ) and \2p 1 / 2 ) functions. Such electroweak radiative corrections for 
the HCI were first obtained in [18] by solving Dirac equation with the Uehling potential 

"-iff-— (^w)^- 

The scheme of calculation of these corrections via PT looks as follows 
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where the notation |2s') and \2p'} means that this functions were obtained as a solutions of the Dirac equation with the 
Uehling potential. The right-hand side of Eq (15) presents PT expansion. We include the second order of PT because of 
its significance. Our results for the matrix elements in the left-hand side of Eq (15) are: (2s| V w \2p) — —4.163986 • 10~ 6 
eV and (2s'\V w \2p') = —4.219566 • 10~ 6 eV. These values are slightly different from the values obtained in [18] because 
of our use of Fermi charge density distribution and more modern data for the nuclear radius. In total the left-hand 
side of Eq (15) is equal to —5, 5580 • 10~ 8 eV. The calculation in the first order of PT in the right-hand side of Eq (15) 
gives —5.5058 • 10~ 8 eV. The evaluation of the second order terms yields the value —5, 2277 • 10~ 10 eV and the sum of 
the first and second orders terms is —5.5581 • 10 -8 eV. Thus the relative discrepancy between the left and right-hand 
sides of Eq (15) is about 2 • 10~ 5 . This is the limit of accuracy of the calculation of the difference between two matrix 
elements (2s|]/ u ,|2p) and (2s'|V^,|2p') that were calculated with accuracy 10~ 7 . 

For our purposes it is interesting to solve the Dirac equation with the PNC weak interaction and to calculate then 
the matrix element of the Uehling potential in a direct way and via PT: 

(2s\V u \2p) =y ^ s Wu\n S )(n S \V w \2p) | y (2 S \V w \np) (np\V u \2p) ^ 

^ E2p — E ns ^ E 2s — Enp 

The results of comparison of the left- and right-hand sides of Eq (16) are given in Table IV. The agreement is about 
1 • 10~ 7 for N w = 1 10. 

It may be interesting also to compare the evaluation of the Uchling-PNC matrix element in a way given by Eq (15) 
and by the application of Eq (16). For this comparison we rearrange the terms in the right-hand side of Eq (16), 
separating out the large contribution, containing small denominator AE^ = E 2s — E 2p . We remind that in our 
calculation we do not include the radiative corrections in the Dirac energy values. Therefore AE 2 ^ only due 
to the use of the non-Coulomb field of the nucleus (Fermi distribution). Using the smallness of the differences 
A£'"' = E ns — E np we replace E ns by E np and vice versa in the right-hand side of Eq (16). Then we arrive at the 
equality 

{2~s\V u \2p) ' (2a|K|2a)] {2s\V w \2p) ^ y {2 S \V u \ns) {n S \V w \2p) | y (2s\V w \np) (np\V u \2p) 

E2s — E2p E2s — Ens E 2p — E np 
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the right-hand side of which coincides with the first two terms in the right-hand side of Eq (15). Evaluation of the 
left-hand side of Eq (17) gives the value —5.5057 • 1CP 8 eV in a good agreement with the first-order PT result in 
Eq (15). 

The second order of PT can be obtained by calculation of the following expression 
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For this expression we obtain the value —5,2239 • 1CT 10 eV. Then the sum of the first and second orders terms is 
—5.5579- 10 -8 eV. Then, the relative accuracy of calculation of the sum of the first and second order terms is the same 
as in the first method: 2 • 1CP 5 . Here we also lose two orders in accuracy due to subtraction of the leading terms in 
Eq (17). In any case we can state emphasize that the accuracy of calculation of any PNC matrix elements is of order 
10~ 7 . Of course, one can lose the accuracy of calculation because of subtraction of two large values. The question of 
the accuracy of calculation of some composite expression should be solved separately in each case. As we understand, 
the calculation of the radiative corrections to the emission PNC amplitude should not contain any subtractions. 

Concluding, we can say that the method developed here presents a powerful tool for the calculation of any corrections 
to the PNC matrix elements in atoms and ions. The recommended scaling parameter is N w = 1 or N w = 10 and the 
relative accuracy for the calculation of PNC matrix elements achieved in all cases is not less than 10 -7 . In particular, 
this approach can be applied to the calculation of the radiative corrections to the emission PNC amplitude with high 
accuracy without any approximations (aZ expansion, Uehling potential approximation and other). Including in our 
scheme the Dirac-Hartree-Fock potential we can calculate these radiative corrections also for neutral atoms, e.g. for 
Cs. It will help to solve finally the problem of the possible deviation of the measured weak charge value Q w from 
that predicted by the Standard Model. 



ACKNOWLEDGMENTS 



The authors thank I.Goidcnko and O.Yu.Andreev for the help with the B-spline codes. They are also grateful to 
G.Soff for drawing their attention to the paper [28]. The financial support by the RFBR grant 02-02-16758 and by 
Minobrazovanje grant E02-3.1-7 is acknowledged. 



[1] CS. Wood et ai, Science 275, 1759 (1997) 

[2] S.C. Bennett and C.E. Wieman, Phys.Rev.Lett. 82, 2484 (1999) 

[3] A. Derevianko, Phy.Rev.Lett. 85, 1618 (2000) 

[4] V.A. Dzuba, C. Harabati, W.R. Johnson and M.S. Safronova, Phys.Rev. A 63, 044103 (2001) 

[5] M.G. Kozlov, S.G. Porsev and I.I. Tupitsyn, Phys.Rev.Lett. 86, 3260 (2001) 

[6] A. Derevianko, Phys.Rev. A 65, 012106 (2001) 

[7] A. A. Vasilyev, I.M. Savukov, M.S. Safronova and HG. Berry, Phys.Rev. A 66, 020101 (2002) 

[8] W.R. Johnson, I. Bednyakov and G. Soff, Phys.Rev.Lett. 17, 233001 (2001) 

[9] A.I. Milstein and O.P. Sushkov, Phys.Rev A 66, 022108 (2002) 

[10] A.I. Milstein, O.P. Sushkov and I.S. Terekhov Phys.Rev.Lett. 89, 28003 (2002) 

[11] M.Yu. Kuchiev, V.V. Flambaum Phys.Rev.Lett. 89, 283002 (2002) 

[12] M.Yu. Kuchiev, V.V. Flambaum J.Phys. B 36, R191-R222 (2003) 

[13] VG. Gorshkov and L.N. Labzowsky, Pis'ma Zh.Eksp.Teor.Fiz. 19, 768 (1974) [JETP Lett. 19, 394 (1974)] 

[14] A. Schafer, G. Soff, P. Indelicate, B. Miiller and W. Greiner, Phys.Rev. A 40, 7362 (1989) 

[15] G. von Oppen, Z.Phys. D 21, 181 (1991) 

[16] V.V. Karasiev, L.N. Labzowsky and A.V. Ncfiodov, Phys.Lett. A 172, 62 (1992) 

[17] R.W. Dunford, Phys.Rev. A 54, 3820 (1996) 



5 



[18] I. Bednyakov, L. Labzowsky, G. Plunien, G. Soff and V. Karasiev, Phys.Rev. A 61, 012103 (1999) 

[19] J. Sapirstein, K. Pachucki, A. Veitia and K.T. Cheng, Phys.Rev. A 67, 052110 (2003) 

[20] L.N. Labzowsky, A. Prosorov, A.V. Shonin, I. Bednyakov, G. Plunien and G. Soff, Annals of Physics 302, 22 (2002) 

[21] P.J. Mohr, G. Plunien and G. Soff, Phys.Rep. 293, 229 (1998) 

[22] S.A. Blundell and N.J. Snyderman, Phys.Rev. A 44, R1427 (1991) 

[23] A. Mitrushenkov, L. Labzowsky, I. Lindgren, H. Persson and S. Salomonson, Phys.Lett. A 200, 51 (1995) 

[24] W.R. Johnson, S.A. Blundell and J. Sapirstein, Phys.Rev. A 37, 307 (1988) 

[25] C. Froese-Fischer and F.A. Parpia, Phys.Lett. A 179, 198 (1993) 

[26] LB. Khriplovich, Parity Nonconservation in Atomic Phenomena, Gordon and Breach, New York, 1991 

[27] D.E. Groom at al, Euro.Phys.J. C 15, 1 (2000) 

[28] A. Schafer, B. Miiller and W. Greiner, Zeit. ffir Phys. A 322, 539 (1985) 

[29] C.W.A. Fletcher, Computational Galerkin Methods, Springer, New York, 1984 

[30] A. Chodos, R.L. Jaffe, K. Johnson, C.B. Thorn and V.W. Weisskopf, Phys.Rev. D 9, 3471 (1974) 

[31] V.A. Yerokhin, A.N. Artemyev, T. Beier, G. Plunien, V.M. Shabaev and G. Soff, Phys. Rev. A 60, 3522 (1999) 



TABLE I. Comparison of the Sommerfeld and spline spectra in absence of V v 
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l s l/2 
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1 9F 1 A 


z6 l/2 


19^7 ^0^8404^0810 
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2.4E-15 
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-1257.395849439792 


-1257.395849439807 


-1.2E-14 


2P3/2 


-1089.611415894742 


-1089.611415894764 


-2.0E-14 


3S1/2 


-539.0933280909740 


-539.0933280909703 
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3fl/2 
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3P3/2 


-489.0370846743426 
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3^5/2 
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4si/ 2 
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-266.3894469008155 


-1.4E-14 


5Si/ 2 


-185.4851885786884 
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TABLE II. Comparison of the exact (point-like nucleus) and spline values for the matrix elements < 2s 1 / 2 |r m |2s 1 / 2 >, 
< 2p 1 / 2 |r m j2p 1/2 >, in absence of V w . 



Matrix element < 2s\r m 2s > in a.u. 
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Spline functions 


Exact functions 


Relative inaccuracy 
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3.428467651648418E-3 


3.428467651648391E-3 


-6.4E-15 
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5.333841373470161E-2 


5.333841373470151E-2 


-2.6E-15 





1.00000000000000 


1.00000000000000 


1.0E-16 


-1 


33.2605595132210 


33.2605595132212 


1.2E-14 


-2 


8415.21963758472 


8415.21963813828 


6.6E-11 




Matrix element < 2p 


\r m \2p > in a.u. 
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Spline functions 


Exact functions 


Relative inaccuracy 
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2.276847838830044E-3 


2.276847838830025E-3 
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-1 
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33.2605595132212 


1.2E-14 


-2 


2542.94240478497 


2542.94246617904 


2.4E-08 



TABLE III. Relative contribution of PNC weak interaction to the energies for the different values of N w (the left-hand side 
of Eq (11) divided by nonperturbative energy value) and the relative contribution of weak interaction to the energies calculated 
by perturbation theory (the right-hand side of Eq (11) divided by nonperturbative energy value). The last column shows the 
relative inaccuracy of spline method from Table I. 



N w 


1 


10 


10 2 


10 3 


10 4 


PT value 


Sspl 


1*1/2 
2#l/2 
2s*i/2 


-2.65322E-13 
-7.13101E-13 
6.36938E-14 


-2.94537E-11 
-3.13915E-11 
6.64008E-12 


-2.95031E-09 
-3.14292E-09 
6.57424E-10 


-2.95046E-07 
-3.14285E-07 
6.57067E-08 


-2.96521E-05 
-3.13759E-05 
6.39177E-06 


-2.95031E-13 
-3.14291E-13 
6.57255E-14 


2.4e-15 
-1.2E-14 
-2.0E-14 



TABLE IV. The matrix elements < 2s"i/ 2 |r m |2pi/ 2 > (m=2, 1,-1,-2) calculated by different methods with the different values 
of the scaling parameter N w = 10 r . The first line in the every box for a certain m and r values corresponds to the direct 
evaluation of the matrix element with B-spline solution of the Dirac equation with PNC weak interaction. The second line 
corresponds to the use of the formula Eq (12) in the text. The PT line corresponds to the perturbation theory evaluation via 
the formula Eq (13) with N w = 1. (in a.u.) The last column presents the results of evaluation of the Uehling matrix element 
with the same methods. 



r 


m 


{2s\V u \2p) eV 




2 


1 


-1 


-2 







-1.971129293262E-8 


-1.858655956533E-7 


-1.914713446350E-7 


-8.871229471602E-2 


1.75928086582129E-6 




-1.971129293278E-8 


-1.858655956541E-7 


-1.914713461356E-7 


-8.871223340185E-2 


1.75928086582128E-6 


1 


-1.971129190577E-7 


-1.858655859887E-6 


-1.914710897826E-6 


-8.871235842834E-1 


1.75928077138616E-5 




-1.971129190573E-7 


-1.858655859884E-6 


-1.914710899371E-6 


-8.871236585343E-1 


1.75928077138618E-5 


2 


-1.971122272187E-6 


-1.858649339848E-5 


-1.914707122846E-5 


-8.87120565454E+0 


1.75927469110662E-4 




-1.971122272187E-6 


-1.858649339848E-5 


-1.914707122268E-5 


-8.87120593112E+0 


1.75927469110664E-4 


3 


-1.970430890847E-5 


-1.857997762638E-4 


-1.914283967935E-4 


-8.86823300844E+1 


1.75866695639621E-3 




-1.970430890846E-5 


-1.857997762637E-4 


-1.914283967449E-4 


-8.86823302428E+1 


1.75866695639620E-3 


4 


-1.904352098651E-4 


-1.795723793234E-3 


-1.874284426839E-3 


-8.58437274571E+2 


1.70059752407138E-2 




-1.904352098651E-4 


-1.795723793234E-3 


-1.874284426869E-3 


-8.58437274521E+2 


1.70059752407139E-2 


PT 


-1.971129278718E-8 


-1.858655943008E-7 


-1.914711189050E-7 


-8.871235465273E-2 


1.75928083527704E-6 
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